# for seekgene CRC human sample
library(Seurat)
library(tidyseurat)
source('00_util_scripts/mod_seurat.R')
library(tidyverse)

# seekgene starsolo velocyto ----------
skg_path <- '~/work/results/seekgene_crc/expr_mat/'

skg_raw <- list.files(skg_path, pattern = 'raw$', include.dirs = TRUE, full.names = TRUE, recursive = TRUE) |>
  str_subset('Gene')

skg_final <- list.files(skg_path, pattern = 'filtered$', include.dirs = TRUE, full.names = TRUE, recursive = TRUE) |>
  str_subset('Gene')

sobj_raw <- skg_raw |>
  set_names(c('crc221125','crc230304','crc230309')) |>
  Read10X() |>
  CreateSeuratObject(min.features = 1)

# save original raw matrix for soupx correction
raw_list <- sobj_raw |>
  SplitObject(split.by = 'orig.ident') |>
  map(GetAssayData, 'count')

sobj_sg <- skg_final |>
  set_names(c('crc221125','crc230304','crc230309')) |>
  Read10X() |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

# drops with cells
celldrop <- sobj_raw |> nest(data = -orig.ident) |>
  pull(data) |>
  map(GetAssayData, .progress = TRUE) |>
  map(DropletUtils::emptyDrops, .progress = TRUE) |>
  map(as_tibble, rownames = '.cell') |>
  map(filter, FDR < .01) |>
  map(pull, .cell) |>
  list_c()

# try plot inflection
sobj_raw <- sobj_raw |>
  CalculateBarcodeInflections(threshold.low = 500, threshold.high = 10000)

sobj_raw |> BarcodeInflectionsPlot()

sobj_raw |>
  pluck('tools', 'CalculateBarcodeInflections', 'inflection_points')

# keep only non-empty drops
sobj_raw <- sobj_raw |>
  filter(nFeature_RNA >= 200 & .cell %in% celldrop)

sobj_raw$mitoRatio <- sobj_raw |> PercentageFeatureSet("^MT-")

sobj_raw <- sobj_raw |>
  filter(mitoRatio < 10)

sobj_raw <- sobj_raw |> quick_process_seurat()

## remove ambient RNA -----
cluster_list <- sobj_raw |>
  as_tibble() |>
  select(.cell, orig.ident, seurat_clusters) |>
  nest(data = -orig.ident) |>
  pull(data) |>
  map(\(x)set_names(x$seurat_clusters, x$.cell))

sobj_list <- sobj_raw |>
  SplitObject(split.by = 'orig.ident') |>
  map(GetAssayData, 'count')

des1 <- adjust_soup_solo(sobj_list[[1]], raw_list[[1]], cluster_list[[1]])

desoup_list <- list(sobj_list, raw_list, cluster_list) |>
  pmap(adjust_soup_solo, .progress = TRUE)

sobj <- desoup_list |>
  map(CreateSeuratObject, min.features = 200, min.cells = 3) |>
  purrr::reduce(merge)

sobj$mitoRatio <- PercentageFeatureSet(sobj, "^MT-")

VlnPlot(sobj, 'mitoRatio')

sobj <- sobj %>%
  filter(mitoRatio < 10)

sobj <- quick_process_seurat(sobj)

## remove doublets ----------
system.time(
sobj.dbl <- sobj |>
  as.SingleCellExperiment() |>
  scDblFinder::scDblFinder(clusters = 'seurat_clusters',
                           samples = 'orig.ident',
                           BPPARAM = BiocParallel::MulticoreParam(tasks = 0, progressbar = TRUE)) |>
  as.Seurat())
  
sobj.dbl |> count(scDblFinder.class)

sobj.dbl |>
  count(seurat_clusters, scDblFinder.class) |>
  ggplot(aes(seurat_clusters, n, fill = scDblFinder.class)) +
  geom_col(position = "fill")

sobj <- sobj.dbl |> filter(scDblFinder.class == 'singlet')

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- mark_cell_type_singler(sobj, ref = hpca, new_label = 'hpca_main')

sobj_raw <- mark_cell_type_singler(sobj_raw, ref = hpca, new_label = 'hpca_main')

DimPlot(sobj, group.by = 'hpca_main',
        cols = DiscretePalette(length(unique(sobj$hpca_main)),
                               palette = 'polychrome'))

sobj_raw |> DimPlot(group.by = 'hpca_main',
                    label = TRUE, label.box = TRUE,
                    cols = DiscretePalette(length(unique(sobj_raw$hpca_main)),
                                           palette = 'polychrome')) 

sobj <- sobj |>
  mutate(genotype = case_when(orig.ident == 'crc221125' ~ 'IT',
                              .default = 'TT')) |>
  filter(!str_detect(hpca_main, 'Endo|Epith|stem'))

sobj_raw <- sobj_raw |>
  mutate(genotype = case_when(orig.ident == 'crc221125' ~ 'IT',
                              .default = 'TT')) |>
  filter(!str_detect(hpca_main, 'Endo|Epith|stem'))

sobj_raw |>
  ggplot(aes(genotype, fill = orig.ident)) +
  geom_bar() +
  coord_flip() +
  labs(fill = 'sample')

write_rds(sobj, 'CRC-I/data/seekgene/crc2305.emtpydrop.rds')

write_rds(sobj_raw, 'CRC-I/data/seekgene/crc2305.emtpydrop.rds')

# guo 2021 sided ---------
guo_mtx <- list.files('CRC-I/data/guo2021sided', pattern = 'mtx', full.names = TRUE)
guo_cell <- list.files('CRC-I/data/guo2021sided', pattern = 'barc',full.names = TRUE)
guo_gene <- list.files('CRC-I/data/guo2021sided', pattern = 'feat',full.names = TRUE)

guo_sparse <- pmap(list(guo_mtx,guo_cell,guo_gene), ReadMtx, .progress = TRUE)

sobj_guo <- guo_sparse |>
  map2(c('left1','left2','left3','right1','right2','right3'), add_name_field) |>
  map(CreateSeuratObject, names.field = 2, min.cells = 3, min.features = 200) |>
  purrr::reduce(merge)

sobj_guo

sobj_guo$mt.ratio <- sobj_guo |>
  PercentageFeatureSet('^MT-')

sobj_guo |> VlnPlot('mt.ratio') 

sobj_guo <- sobj_guo |> filter(mt.ratio < 10)

sobj_guo <- sobj_guo |>
  quick_process_seurat()

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj_guo <- mark_cell_type_singler(sobj_guo, hpca, new_label = 'hpca_main')

sobj_guo |> DimPlot(group.by = 'hpca_main',
                    cols = DiscretePalette(11, palette = 'polychrome'))

immune_guo <- sobj_guo |>
  filter(str_detect(hpca_main, 'B|DC|Macro|Mono|Neutro|NK|T'))

immune_guo |> DimPlot(group.by = 'hpca_main',
                    cols = DiscretePalette(7, palette = 'polychrome'))

immune_guo |> VlnPlot('ACTB', group.by = 'hpca_main')

immune_guo |>
  mutate(genotype = ifelse(str_detect(orig.ident, 'left'), 'IT', 'II')) |>
  write_rds('CRC-I/data/guo2021sided/guo2021_immune.rds')

immune_guo <- read_rds('CRC-I/data/guo2021sided/guo2021_immune.rds')
immune_guo |> count(genotype)

immune_guo <- immune_guo |>
  quick_process_seurat(pcs = 25, res = 2)

immune_guo <- immune_guo |>
  mark_cell_type_singler(ref = sce_10x,
                         sc_ref = TRUE,
                         new_label = 'zhang2020_main')

immune_guo |> DimPlot(group.by = 'zhang2020_main') +
  ggtitle('Major tumor immune cell groups')

# li2021liver ------------
sobj_li <- Read10X('CRC-I/data/li_2021_liver')

sobj_liv <- sobj_li |>
  CreateSeuratObject(names.field = 2, min.cells = 3, min.features = 200)

sobj_liv |> dplyr::count(orig.ident)

sobj_liv <- sobj_liv |> mutate(tissue = str_remove(.cell, '.+_'))

sobj_liv$mitoRatio <- sobj_liv |> PercentageFeatureSet("^MT-")

sobj_liv |> VlnPlot('mitoRatio')

sobj_liv <- filter(sobj_liv, str_detect(tissue, 'CRC') & mitoRatio < 10)

meta_liv <- "COL07	IT
COL12	II
COL15	II
COL16	IT
COL17	IT
COL18	II
" |> read_delim(col_names = c('orig.ident','genotype'))

sobj_liv <- sobj_liv |> left_join(meta_liv)

sobj_liv |> count(genotype)

sobj_liv <- sobj_liv |> quick_process_seurat()

sobj_liv <- sobj_liv |>
  FindVariableFeatures() |>
  ScaleData() |>
  RunPCA() |>
  RunHarmony('orig.ident') |>
  RunUMAP(reduction = "harmony", dims = 1:20) |>
  FindNeighbors(reduction = "harmony", dims = 1:20)

sobj_liv |> write_rds('CRC-I/data/li_2021_liver/li2021crc.rds')

sobj_liv <- read_rds('CRC-I/data/li_2021_liver/li2021crc.rds')

system.time(sobj_liv <- sobj_liv |>
              FindClusters(algorithm = 4,
                           resolution = 1,
                           method = 'igraph'))

sobj_liv |> DimPlot()

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj_liv <- mark_cell_type_singler(sobj_liv,
                                   hpca,
                                   new_label = 'hpca_main')

sobj_liv |> DimPlot(group.by = 'hpca_main',
                    cols = DiscretePalette(10))

sobj_liv <- sobj_liv |> filter(!str_detect(hpca_main, 'Endo|Epi|stem'))

sobj_liv |> write_rds('CRC-I/data/li_2021_liver/li2021crc.rds')

sobj_liv <- read_rds('CRC-I/data/li_2021_liver/li2021crc.rds')

sobj_liv <- quick_process_seurat(sobj_liv)


